Mathematical Biosciences

Cancer-associated fibroblasts (CAFs) are a central component of the tumor microenvironment that facilitate a supportive niche for cancer progression and metastasis. Experimental evidence suggests that CAFs can facilitate estrogen-independent tumor growth, thereby reducing the efficacy of anti-hormonal therapies. Understanding and quantifying the complex interactions between tumor cells, hormonal signalling, and the microenvironment are crucial for designing more effective and individualized treatment strategies. We propose a mathematical framework to explore the influence of CAFs on ER+ breast cancer progression and to evaluate strategies to mitigate their impact. We develop a system of nonlinear ordinary differential equations that substantiates the experimental observations by providing a mechanistic basis for the role of CAFs in regulating estrogen-independent growth dynamics. We then employ optimal control theory to evaluate distinct therapeutic approaches involving monotherapy or combinations of: (i) inhibition of tumor-to-CAF signaling, (ii) inhibition of CAF-to-tumor proliferative signaling, and (iii) endocrine therapy. Taken together, our results demonstrate that CAF-targeted strategies can enhance treatment efficacy across various estrogen dosing regimens. Our study provides new insights into the potential of CAF as a therapeutic target that could help to improve existing approaches for endocrine therapies.

We present an approach to analysing cell homeostasis using a bond graph modelling approach that ensures that the conservation laws of physics (conservation of mass, charge, and energy, respectively) are satisfied for the interdependent biochemical, electrical, mechanical, and thermal energy storage mechanisms operating within the cell. We apply the bond graph approach to several cell membrane transport mechanisms and then consider how physics constrains intracellular electrolyte homeostasis for enterocytes (the epithelial absorptive cells of the gut). The model includes the electrogenic sodium-potassium ATPase pump (NKA), the glucose transporter (GLUT2), and an inwardly rectifying potassium channel, all in the basolateral membrane, and the electrogenic sodium-driven glucose transporter (SGLT1) in the apical membrane. Glycolysis converts the imported glucose to ATP to drive NKA. For specified levels of sodium, potassium, and glucose in the blood, the model demonstrates how enterocytes absorb sodium and glucose from the gut and transfer glucose to the blood while maintaining the membrane potential and homeostasis of intracellular sodium and potassium. The Gibbs free energy available from the ATP hydrolysis ensures that the cell operates as a sodium battery with a high external to internal ratio of sodium concentration in order to provide the energy for many other cellular transport processes. We show that the 3:2 stoichiometry of Na+/K+ exchange in NKA, coupled with 2:1 Na+/glucose cotransport in SGLT1, a 1:2:2 ratio between glucose consumption and ATP and water production in glycolysis, and K+ and glucose efflux through Kir and GLUT2, respectively, provides a balanced system that maintains homeostasis of intracellular Na+, K+, glucose, ATP and water, and homeostasis of the membrane potential, under varying levels of transport of glucose from the gut to the blood. We also show how the flux expressions for SLC transporters, ATPase pumps and ion channels can all be expressed in a consistent and thermodynamically valid way.

Vaccine-acquired immunity plays an important role in controlling the spread of many infectious diseases; however, vaccine efficacy can diminish over time. This work uses a mathematical model to study the effects of waning vaccination-acquired immunity on infection incidence. With an SEIR-type compartmental model that considers both vaccinated and unvaccinated populations (and their mixing), we present mathematical conditions under which vaccinated individuals drive ongoing growth in infections, i.e., over half of the daily incidence arises from vaccinated individuals. Analysis of a mathematical model of COVID-19 spread in the state of Colorado suggests how and for what duration vaccinated individuals could have sustained such growth. Importantly, our model demonstrates that, despite potential for brief vaccinated-driven periods of growth in infections, which occur among unvaccinated-driven periods of growth in infections, increased vaccination coverage always reduces total cases and total hospitalizations. This work provides insight into how waning immunity in vaccinated populations can contribute to ongoing infection incidence and demonstrates the value of complementary interventions to prevent disease spread in vaccinated populations.

Following the suggestion of L. F. Jaffe [1] we propose a mathematical model of fast calcium induced calcium influx waves (CICI Waves). They can propagate at relatively high speeds (up to 1300 micrometers/s). According to [1], they propagate due to a mechanochemical interaction of actomyosin network with the cell membrane. The local stretching of the membrane caused by actin filaments opens mechanically operated ion channels resulting in the influx of calcium to the cell. Moreover, stretching a cells membrane at one point opens nearby stretch activated calcium channels because the mechanical force is relayed by the actin filaments interconnected by myosin bridges. The number of bridges as well as filament density increases with calcium concentration, causing the contraction of the actomyosin network. Thus, the force acting on the membrane from tangled actin filaments is transmitted ahead of the moving front of the calcium concentration. As a result, the ion channels are opened even before the signal of calcium reaches them. This leads to much larger propagation speed of CICI waves in comparison with calcium induced calcium released (CICR) waves, where the wave is sustained by the diffusion of calcium and autocatalytic release of calcium from the internal stores (e.g. endoplasmic reticula).

The placenta is a fundamental organ for human reproduction, facilitating fetal growth via the exchange of oxygen, nutrients and waste products between mother and fetus, by means of a dense network of fetal villi bathed in maternal blood that flows through the intervillous space (IVS). However, despite its role in adverse pregnancy outcomes associated with impaired maternal-fetal transfer, the influence of placental structure on maternal haemodynamics and the delivery of oxygen and nutrients to the developing fetus is not well understood. This study employs computational fluid dynamics within physiologically informed 2D and 3D whole-organ-scale placental geometries, whose features are informed by recent ex vivo experimental and {micro}CT data, to examine comprehensively the influence of placental anatomy on maternal flow and transport in the IVS. In particular, we consider in detail the impact of the number and placement of maternal decidual arteries and veins that supply and drain the IVS, the location and height of so-called septal walls that loosely separate the placenta into functional units (cotyledons) and sub-units (lobules), and the density of the fetal villous trees as reflected in the rate of uptake of dissolved solutes from the maternal blood and the resistance to flow. We first exploit the computational efficiency of simulation in a representative 2D geometry to study in detail the sensitivity of haemodynamic markers to these parameters. These results guide our 3D study which reveals that the flow, transport and oxygen uptake are strongly influenced by placental structure, and exposes the vein-to-artery ratio as a key indicator of placental efficiency, regulating a trade-off between a preferential maternal flow environment and fetal oxygen uptake. Conversely, the location and height of the septal walls, a feature that is not well studied, have minimal systematic impact on the macroscopic haemodynamic and transport measures considered here. We also introduce a reduced model, for which analytical progress can be made, and demonstrate its utility in exposing key drivers of maternal-fetal transport. Author summaryIn our study, we explored how placental structure affects maternal blood flow and the delivery of oxygen to the fetus. The placenta is a complex and vital organ, but we do not fully understand how its morphology influences its functional efficiency. This is a critical gap in our knowledge as placental blood flow is implicated in serious pregnancy complications such as fetal growth restriction and pre-eclampsia. To investigate the effect of structure, we used detailed 2D and 3D models of the placenta. These models allowed us to simulate maternal blood flow and oxygen transport therein while changing various structural features, such as the number and location of placental veins and the density of the intervillous space. We found that the ratio of veins to arteries is a key factor in determining the maternal blood flow speed and the amount of oxygen the fetus receives. A higher number of veins relative to arteries helps create the slow-flow environment necessary for a healthy pregnancy, but an excessive ratio can reduce the efficiency of oxygen uptake. Our findings suggest that the vein-to-artery ratio could be an important indicator of placental health, offering a new perspective for understanding and enabling effective early intervention treatments of pregnancy complications.

This study presents a theoretical and mathematical framework for understanding the dynamical behavior of infectious disease spread using a compartmental modeling approach. The proposed model incorporates memory effects to capture temporal dependencies that are not adequately represented by classical formulations. Qualitative analysis is employed to investigate the stability properties of the system and the role of key mechanisms in shaping long term dynamics. Publicly available surveillance information is used only to illustrate the consistency of the model behavior with observed trends. The results highlight the value of memory based modeling structures for describing complex biological processes and provide a general mathematical perspective for studying epidemic dynamics.

Evolutionary therapies regulate heterogeneous populations by altering selective pressures through treatment sequences in cancer and infections. This letter develops an invariant-set framework for treatment-induced containment based on positive triangular invariant sets. For periodically switched systems, sufficient conditions are derived for the existence of such invariant regions. Robustness with respect to mutation is established by showing that the invariant simplex persists under small perturbations of the subsystem matrices. In the two-phenotype case, the analysis yields an explicit mutation threshold that separates regimes in which therapy cycling maintains containment from regimes in which mutation can enable evolutionary escape. Simulations illustrate the geometry of the invariant sets and the role of mutation and dwell time in containment robustness.

Human immunodeficiency virus (HIV) persistence remains a major barrier to cure due to the existence of long-lived latent reservoirs that evade immune clearance and persist despite combination antiretroviral therapy (ART). Although ART effectively suppresses viral replication, treatment interruption often leads to rapid viral rebound originating from these latent reservoirs. In this study, we develop a deterministic mathematical model describing the in vivo dynamics of HIV infection incorporating uninfected CD4+ T cells, infected cells, latent reservoirs, deep latent reservoirs, and infectious and non-infectious virions, while explicitly accounting for the therapeutic effects of reverse transcriptase inhibitors (RTIs), protease inhibitors (PIs), and Tat transcription inhibitors. Analytical results establish positivity and boundedness of solutions and derive the effective reproduction number Re using the next-generation matrix approach. Stability analysis shows that the virus-free equilibrium is locally asymptotically stable when Re < 1, while viral persistence occurs when Re > 1. Numerical simulations were performed to investigate therapy interactions, viral rebound following treatment interruption, and the impact of drug efficacy on viral set-points and latent reservoir dynamics. To further explore therapy interactions, three-dimensional viral set-point surfaces and heat maps were generated to examine how combinations of infection inhibition, viral production inhibition, and transcriptional inhibition influence viral dynamics. The simulations reveal that Tat inhibition suppresses viral transcription, thereby reducing the transition of infected cells into productive infection and limiting viral propagation when combined with conventional ART mechanisms. The therapy parameter planes further demonstrate that strong transcriptional inhibition promotes the transition of infected cells into deep latency, supporting the emerging block-and-lock strategy for functional HIV cure. In addition, a three-dimensional eradication boundary surface and therapy cube were constructed to identify regions of parameter space where Re < 1, corresponding to successful viral control. These visualizations show that viral eradication is unlikely when therapies act independently but becomes achievable when multiple therapeutic mechanisms act simultaneously. Overall, the results highlight the critical role of transcriptional inhibition through Tat-targeting therapies in complementing existing ART regimens. By simultaneously suppressing viral replication and promoting deep latency, Tat-based combination strategies may significantly reduce viral rebound and contribute to long-term functional control of HIV infection.

Atherosclerotic plaques are fatty deposits in arterial walls and a major cause of heart attacks and strokes. Macrophage proliferation triggers plaque growth and instability, but the exact conditions that cause stable plaques to become unstable remain unclear. To provide an insight into the conditions for this transition, we apply bifurcation analysis to the lipid-structured atherosclerosis model proposed by Chambers et al. (Bull Math Biol 86(8):104, 2024). Our main contribution is that the reduced dynamics of the system remain meaningful even beyond previously identified limits of validity. Furthermore, along with numerical bifurcation methods, the use of fast-slow analysis, combined with Fenichels theory, identifies a saddle-node bifurcation at infinity. A sharp threshold exists where macrophage proliferation balances emigration. Below this balance, the system stabilises in a biologically reasonable state; contrary to above it, macrophage numbers and lipid load grow unboundedly, triggering instability and runaway inflammation. Trends in determinant and eigenvalues also support this threshold. Parameter scans and heatmaps demonstrate that increased proliferation or reduced emigration enhances the number of macrophages and the lipid content of the necrotic core. Efferocytosis rate modulates downstream severity but does not shift the primary threshold. These findings reconcile conflicting results on macrophage proliferation, demonstrating that it is protective when emigration sufficiently balances this process. In other words, co-targeting reduced macrophage proliferation and enhanced emigration could help maintain plaque stability and reduce the risk of acute cardiovascular events. While this remains a theoretical recommendation, it offers a potential therapeutic strategy that authorises further investigation in experimental and clinical settings.

Intratumor heterogeneity (ITH) arises from the combined effects of genetic alterations, clonal interactions, and environmental constraints, and plays a central role in therapeutic resistance and disease progression. While ITH has been extensively documented in empirical tumor data, the scientific debate regarding the biological mechanisms underlying this heterogeneity remains complex, highlighting the need for cancer evolution models that are sufficiently flexible and sophisticated to reproduce the observed behaviors and to give insights on the unobserved ones. Here, we present a stochastic modelling framework for tumor evolution that integrates genotypic inheritance with phenotype driven functional traits and resource mediated competition. Mutational events are associated with functional capabilities such as altered proliferation, increased mutation rates, limit evasion potential or enhanced control over shared resources, allowing multiple genotypes to converge on similar phenotypes. The model explicitly tracks subclonal lineages while incorporating environmental constraints that modulate growth and competition.The framework is defined through a mathematically rigorous construction and is accompanied by an efficient simulation algorithm. To facilitate exploration and reproducibility, we provide an open-source graphical user interface that allows users to configure model parameters, run simulations, and inspect clonal genealogies and population dynamics without requiring direct interaction with the underlying code. Using this model, we illustrate how ecological feedbacks can shape clonal dynamics over time, supporting an interpretation in which early tumor growth is dominated by stochastic expansion, while later evolution increasingly reflects selection for traits that alleviate environmental constraints. Rather than constituting a new evolutionary paradigm, this behaviour demonstrates how well-documented biological patterns can emerge naturally from a unified stochastic and ecological description. Overall, our approach offers a flexible and extensible platform for investigating how chance, functional traits, and environmental interactions jointly govern tumor heterogeneity. Author summaryNot all cancerous cells are created equal: inside the same tumor, different populations of cells exist at the same time, fighting for the same resources and influencing the way the disease evolves and reacts to treatments. These groups of cells have different behaviour and abilities thanks to different genetic mutations, which might give them an advantage or bring their population to disappearance. We have built a mathematical model that mimics the evolution of a tumor over time, simulating a competition between its different populations of cells. Our simulated experiments show that tumors evolve in two distinct phases: at first, cells that grow and divide more quickly have an advantage. Once the space and nutrients are limited, cells that can survive with fewer resources have an advantage and can potentially take over the race. We use these simulations to argue that the evolution of a tumor doesnt depend on the shape of the space it expands in, but rather on the availability of nutrients.

Visceral leishmaniasis (VL) is considerably more severe among individuals infected with human immunodeficiency virus (HIV), leading to higher parasite loads, frequent relapse, and increased mortality. To examine the epidemiological interaction between the two diseases, we develop a comprehensive VL-HIV co-infection model that incorporates transmission pathways, treatment effects, and relapse dynamics. The model is parameterized using real-time data from Bihar, India, including monthly VL-only and VL-HIV co-infected cases and annual HIV prevalence data. Our analysis shows that HIV infection drives the resurgence and persistence of VL even in settings where VL alone would not sustain transmission, underscoring the amplifying effect of HIV-induced immunosuppression on VL dynamics. We further demonstrate that increasing HIV treatment coverage substantially reduces co-infection prevalence and lowers VL relapse rates. Numerical simulations and optimal control analysis highlight the effectiveness of integrated intervention strategies that combine awareness, treatment enhancement, and vector control. Overall, this study emphasizes the need for coordinated VL and HIV control programs and provides data-driven guidance for designing sustainable intervention strategies in endemic regions.

Variability in gene expression among single cells and growing cell populations can arise from the stochastic nature of protein synthesis, which often occurs in random bursts. This study investigates the variability in the expression of a growth-sustaining protein, whose concentration is regulated by a negative feedback loop due to cell growth-induced dilution. We model the distribution of protein concentration using a Chapman-Kolmogorov equation for single cells and a population balance equation for growing cell populations. For single cells, we derive an explicit solution for the protein concentration distribution in state space and represent it as a Bessel function in Laplace space. For growing populations, we find that the distribution satisfies a Heun differential equation with singular boundary conditions. By addressing the central connection problem for the Heun equation, we quantify the population-level protein distribution and determine the Mathusian parameter, which characterizes population growth. This work provides a comprehensive analytical framework for understanding how stochastic protein synthesis impacts gene expression variability and population dynamics.

Nipah virus (NiV) is a sporadic yet extremely deadly zoonotic pathogen, with reported case fatality rates of 40%-75% in impacted areas. Prolonged incubation, documented relapse, and delayed-onset encephalitis following apparent recovery indicate that NiV dynamics are influenced by intricate temporal processes. However, mechanistic contributions of these processes to epidemic persistence remain poorly understood. In this study, we develop and analyze a delay differential equation model for NiV transmission that explicitly incorporates incubation delay, relapse, and post-recovery delay effects. We compute a primary-transmission reproduction threshold (R0), characterize the disease-free and endemic equilibria, and analyze their stability, including delay-induced Hopf bifurcations. We show that relapse modifies the endemic-equilibrium existence condition, so an endemic equilibrium is not determined solely by the classical threshold criterion R0 = 1. We calibrate the model to NiV incidence data from Bangladesh (2001-2024) and perform simulations and sensitivity analyses to evaluate the effects of relapse and delays across epidemiological scenarios. Results indicate that sustained oscillations occur only under hypothetical parameter regimes, suggesting that delay-induced periodic outbreaks are unlikely under empirically informed conditions. Scenario analyses demonstrate that relapse and encephalitis-related delays predominantly influence post-peak dynamics, while incubation delay alters the time and intensity of the epidemic peak. We also introduce a relapse-driven replenishment fraction to quantify contribution of relapse to continued transmission, demonstrating its growing significance following the first outbreak peak. Overall, our results identify relapse as a key mechanism for epidemic persistence and underscore the importance of incorporating relapse and biological time delays into epidemiological modeling and public health strategies.

Increasing emissions of CO2 into the atmosphere are causing ocean acidification, threatening calcifying organisms. In this study, we model the physiological responses of coccolithophorids to acidification to understand the ecological and evolutionary outcomes of a system in interaction with zooplankton. Assuming a trade-off between growth and protection against grazing, we show that calcification has bivalent effects on transfers between two trophic levels and that acidity can strongly alter energy transfers. Taking into account the evolution of calcifying phenotypes in response to acidification, we show that the system outcome contrasts with previous results. While the effect of evolution depends on how calcification affects grazing, it nevertheless follows that acidification leads to a decrease in calcifying capacity. This evolutionary decrease may be progressive, but can also lead to tipping points where abrupt shifts may occur. Such a counter-selection of calcification in turn affects ecosystem functioning, enhancing energy transfers within the system and modifying carbon fluxes. We discuss how such eco-evolutionary changes may impact food webs integrity, carbon sequestration into the deep ocean and therefore endanger the carbon pump stability.

Intra-tumor heterogeneity in proliferation rates fundamentally influences cancer progression and treatment resistance. To investigate how continuous phenotypic variation shapes eco-evolutionary dynamics, we develop a phenotype-structured partial differential equation framework that explicitly models proliferation het-erogeneity as a dynamic trait distribution. Our model integrates three key biological principles: (1) phenotypic diffusion capturing heritable variation in proliferation rates, (2) global resource competition enforcing density-dependent growth constraints, and (3) an experimentally grounded life-history trade-off linking elevated proliferation to increased mortality. Using adaptive dynamics, we derive the optimum proliferation rate in a growing tumor, showing that the optimal phenotype dynamically shifts toward slower proliferation as tumors approach carrying capacity. We perform in silico treatment simulations for four different treatment regimes (pan-proliferation, low-, mid-, and high-proliferation targeting) to show how therapeutic selective pressures reshape fitness landscapes. While all treatments slow down tumor growth, they induce divergent evolutionary trajectories: low- and mid-proliferation targeting enrich fast-proliferating clones, whereas high-proliferation targeting selects for slower phenotypes. We connect these dynamics with changes in mean proliferation rates during and after treatment. We use adaptive dynamics to explain the shifts in mean proliferation rate during treatment, showing how each regimen alters the maximum fitness proliferation rate. Our work establishes a predictive, evolutionarily grounded framework for understanding how therapy reshapes tumor proliferation landscapes, offering a mechanistic basis for designing strategies that anticipate and counteract adaptive resistance.

We use simple mathematical models to explore the factors that influence the evolutionary emergence of mpox to a pathogen capable of sustained human to human transmission that poses a global threat. Smallpox eradication followed by the discontinuation of immunization with vaccinia has led to a decline in the level of population immunity against related poxviruses such as mpox. This decline in immunity results in an increase in both the number of spillovers and the extent of human to human transmission. We find that increases in transmissibility of mpox between humans have a much greater effect on the probability of evolutionary emergence compared with increases in the number of zoonotic spillovers. We suggest that while mpox only needs to have a reproductive number slightly greater than one to become endemic, subsequent adaptation is likely to further increase its transmissibility in the human population. As a consequence a much higher level of vaccination (or other intervention) is needed to control the pathogen after its evolutionary emergence compared with what is needed to prevent it from emerging in the first place.

Tetanus is a severe disease of the nervous system, transmitted through bacteria in the environment. In the absence of medical attention, case fatality rates are extremely high. Despite progress towards maternal and neonatal tetanus elimination targets, tetanus remains a serious public health problem. Routine infant and maternal vaccination have contributed to considerable reduction in cases and deaths from tetanus globally. However, protective immunity wanes over time. To increase duration of protection, the World Health Organization recommends three diphtheria-tetanus-pertussis-containing vaccine booster doses be given in early childhood, childhood, and adolescence. Evidence to support country-level decision-making about the introduction of these booster doses is critical. We have developed a novel age-structured, deterministic compartmental model of tetanus transmission and vaccination. The model is driven by environmental transmission and incorporates interventions like hygiene and safe birth practices to reduce the magnitude of environmental transmission. It explicitly models vaccination, separating each dose of the primary series, booster doses, and maternal vaccination to capture dose-specific effectiveness and duration of protection. The model captures heterogeneous immunity profiles by dose and age, and the cumulative nature of vaccine-derived protection. The immune dynamics follow the patterns described in literature and can replicate seroprevalence studies, although the exact characterisation of immunity in the literature still has gaps. This model presents a substantial advancement on previously published models and is well positioned to inform tailored vaccination strategies to reduce neonatal and non-neonatal tetanus.

The transbilayer distribution of plasma membrane cholesterol remains uncertain despite repeated analysis. We propose a new mechanism driving cholesterol sidedness: sterols form simple stoichiometric associations with phospholipids. Our model postulates that the phospholipids in the plasma membrane bilayer are fully complexed with cholesterol. The cholesterol in each leaflet is then the product of the abundance of its phospholipid and its sterol stoichiometry. Notably, lipid affinities are not relevant. Applying literature values for the composition, abundance and sterol stoichiometry of the phospholipid in each leaflet, the model predicts that two-thirds of the cholesterol in the human erythrocyte membrane bilayer is located in its outer leaflet, an exofacial to endofacial ratio of 2:1. The model also predicts that the overall cholesterol content of the bilayer is [~]0.75 mole/mole phospholipid, in agreement with literature values. Furthermore, our analysis suggests that the areas of the two membrane leaflets are about the same. The concordance of prediction with observation validates the model and the values used for the parameters. The sterol in the exofacial leaflet of the plasma membrane of any cell is predicted to exceed that on its contralateral side when its phospholipids have a higher sterol stoichiometry and are fully complexed. SynopsisWe propose that the transbilayer distribution of cholesterol in the plasma membrane bilayer is determined by its complexation with the phospholipids in the two leaflets. Because the complexes are homeostatically filled to stoichiometric equivalence, leaflet cholesterol is given by the abundance of its phospholipids multiplied by its sterol stoichiometry. The model predicts that two-thirds of the cholesterol in the human erythrocyte membrane bilayer resides in the outer leaflet. It also predicts the cholesterol content of the bilayer as a whole.

Calibration of mechanistic cardiovascular models is a central barrier to their use in population analysis and patient-specific simulation, particularly in settings where key physiological variables are unobservable and multiple parameter combinations can reproduce the same haemodynamic targets. In this work, we present Embedded Gradient Descent (EGD), a calibration framework for ODE-based lumped-parameter cardiovascular models in which selected physiological parameters are promoted to dynamic states and driven toward prescribed targets through embedded controller equations. By exploiting the qualitative structure of the governing equations, EGD enforces physiologically consistent parameter-variable relationships, yielding unique calibrated solutions that are robust to initial conditions and scale efficiently with model complexity. The framework is demonstrated using a mechanistic cardiovascular model to generate virtual paediatric populations spanning normal physiology and two septic shock phenotypes (warm and cold shock), achieving low residual error across pressures, flows, and compartmental volumes. The resulting parameter distributions are consistent with known haemodynamic adaptations in paediatric sepsis, including alterations in vascular resistance, compliance, cardiac elastance, and effective blood volume. Importantly, persistent calibration residuals arise only when target combinations are structurally incompatible with the model, providing an explicit and interpretable diagnostic of feasibility limits rather than an optimisation failure. These results establish EGD as a general, scalable calibration strategy for mechanistic cardiovascular models and a practical foundation for virtual population generation and future patient-specific digital twin applications in critical care. NEW & NOTEWORTHYThis study introduces a novel, embedded gradient descent calibration framework that enables scalable generation of mechanistically interpretable virtual populations of patients from ODE-based cardiovascular models. By treating parameter inference as a dynamical extension of the governing equations and calibrating directly against cycle-derived physiological targets, the method preserves physiologically meaningful parameter-variable relationships. Applied to paediatric sepsis, the framework reproduces warm and cold shock phenotypes while exposing infeasible target combinations, while providing efficient calibration and physiological insight.

Rhythmic mood fluctuations have long been linked to circadian ([~]24 h) timing, but how this physiological rhythmicity relates to pathological affective cycling such as bipolar episode recurrence remains unclear. Here we introduce a slow-fast dynamical model with multistability, motivated by hypothalamic-pituitary-adrenal (HPA) axis feedback. We build on a reduced slow-fast formulation with a slow endocrine variable and a fast affective variable, and reshape the fast nullcline to yield four stable fixed points and allow the model to distinguish between normal diurnal mood variation and pathological, depression-like/mania-like extremes. A sinusoidal circadian drive promotes regular alternation within the physiological pair, while temporally correlated fluctuations modeled as an Ornstein-Uhlenbeck process trigger probabilistic escapes. Simulations show that weakened circadian amplitude increases the probability of transitions into pathological attractors and produces prolonged dwell times in affective extremes. Small geometric biases in the nullcline can generate predominant polarity toward depressive or manic episodes. The model provides a conceptual framework linking circadian forcing, persistent stochastic perturbations, and multistability, and generates testable predictions for how circadian disruption destabilizes affective trajectories in mood disorders.